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ABSTRACT 


This  work  was  motivated  by  the  problem  of  obtaining  a  smooth  density  func¬ 
tion  over  a  geographical  region  from  data  aggregated  over  irregular  subregions. 

*  t 

Minimization  of  a  family  of  roughness  criteria  given  ’’volume’’  data  lead  to  smooth 
multivariate  functions  -  Laplacian  histospl ines ,  having  a  certain  order  of  the 


iterated  Laplacian  of  constant  value  in  each  of  the  subregions  and  satisfying 


natural  boundary  conditions  on  the  boundary  of  the  region.  For  inexact  data,  e.g 


in  case  of  estimating  an  underlying  density  given  counts  of  events  by  subregions, 


laplacian  smoothing  histospl ines  are  constructed , analogous  to  smoothing  splines 


in  the  univariate  case,  and  a  method  for  choosing  the  smoothing  parameter  is 


For  both  cases  of  exact  and  inexact  data,  modified  roughness  criteria,  inde 


pendent  of  the  region,  are  discussed,  and  results  known  for  point-evaluation  data 


are  extended  to  the  case  of  aggregated  data 


AMS  (MOS)  Subject  Classifications  41A63,  41A15 


Key  Words:  Hi stospl ines ,  Laplacian  histosplines,  Volume  matching  surfaces 


Bounded  domains.  Smoothing  histosplines.  Elliptic  boundary  value 


problems,  Iterated  Laplacian,  Aggregated  data 


Work  Unit  Number  6  -  Spline  Functions  and  Approximation  Theory 


On  sabbatical  from  Department  of  Mathematical  Science&^Tel-Aviv  University,  Israel 


SI  :NiriCANCK  AND  {EXPLANATION 


•  consider  the  problem  of  obtaining  a  smooth  density  function  when  only  aqqre- 
!4te-t  data  is  available.  For  example,  suppose  that  population  census  is  given  by 
.  ur-aucratic  region  (say,  state)  and  it  is  desired  to  obtain  a  smooth  function  f (x,y) 
intended  to  be  an  estimate  of  the  population  density  at  location  (x,y).  We  obtain 
the  "smoothest"  f  such  that  the  volume  of  f  over  each  region  coincides  with  the 
known  population  size  in  that  region.  Our  measure  of  the  roughness  of  f  is 
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where  the  integral  is  taken  over  the  region  of  interest.  Some  other  rouqhness 
measures  are  also  considered.  The  solution  f  is  characterized  as  the  solution 
to  a  certain  boundary  value  problem.  We  then  modify  the  roughness  criteria  by  taking 
the  integral  over  the  infinite  plane.  The  solution  to  the  modified  problem  can  be 
displayed  explicitly  and  a  computable  approximate  solution  is  obtained. 

We  also  solve  the  problem  of  obtaining  a  smooth  density  when  the  data  must  be 
considered  to  be  inexact,  for  example,  when  it  is  count  data  for  some  rare  disease. 

In  this  case  one  usually  does  not  want  the  volumes  of  f  over  each  region  to 
match  the  data  exactly  but  to  be  near  it.  There  is  a  parameter  controlling  a 
tradeoff  between  the  smoothness  of  f  and  its  deviation  from  the  data,  and  we  show 
how  to  choose  it. 

We  hop>e  that  these  results  provide  a  first  step  in  the  development  of  method; 
for  the  construction  of  surfaces  from  aggregated  data. 
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ON  THE  ESTIMATION  OF  FUNCTIONS 
OK  SEVERAL  VARIABLES  FROM  AGGREGATED  DATA 


Nira  Dyn  and  Grace  Wahba 

1.  Introduction 

The  work  in  this  paper  is  motivated  by  the  following  problem:  Incidence 
rates  of  certain  types  of  cancer  are  known  to  vary  geographically,  for  example, 
persons  living  in  areas  with  higher  exposure  to  sunshine  are  more  likely  to  get 
skin  cancer  than  those  in  more  northerly  regions.  Data  on  population  density 
and  disease  occurrence  is  typically  collected  by  bureaucratic  subdivision.  It 
is  desired,  from  this  aggregate  data,  to  obtain  an  estimate  pCx^x^)  of  the 
probability  pCx^.x.,)  that  a  person  living  at  (x^.Xj)  will  contract  the 
disease  in  a  given  year.  Contour  map  representations  of  p  can  then  be  used 
to  visually  look  for  geographic  patterns  in  p,  and  for  apparent  correlations 
with  other  geographically  varying  variables. 

For  concreteness,  we  consider  data  reported  by  state.  Let  Cl  represent 
the  contiguous  48  states  of  the  U.S.,  and  CL  the  ith  state.  If  uCx^.x^)  is 
the  population  density  at  point  (x^.x^)  (we  pretend  this  is  well  defined),  then 
the  expected  number  of  cases  of  our  subject  disease  in  state  i  is 

u i  “  /  p(x1,x2)  u(Xj,x2)dXjdx2  . 

“i 

The  population  s.  *  j  u(x  ,x  )dx.dx  of  state  i  is  assumed  to  be  known 

Cli  12  12 

exactly.  The  population  of  further  subdivisions,  e.g.,  counties,  can  also  be 
assumed  to  be  known  exactly.  In  a  particular  year  the  number  Z^  of  cases 
actually  occurring  in  Cl^  is  reported.  If  p  is  very  small,  then  Z^  may  be 
modelled  as  a  Poisson  random  variable  with  mean  From  this  data  it  is  desired 

• 
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to  estimate  ptx^.x^),  (Xj.x,)  <.  ft.  We  will  do  this  by  first  estimating  u(Xj.x,) 
using  only  the  population  data  (s^),  and  then  estimating  g(Xj,x.,)  H  p(x^  .x^)  u(x  j  ,x., 
using  the  disease  count  data  {7^}.  The  estimate  of  p  is  then  the  quotient  of 
these  two  estimates.  For  notational  convenience  we  suppose  that  population  data  is 
aggregated  at  the  same  level  (i.e.  state)  as  the  disease  count  data. 

It  is  possible  to  obtain  heuristically  reasonable  estimates  of  u  and  g  by 
assuminq  that  they  are  "smooth"  in  some  sense,  namely  by  minimizing  certain  measures 
of  roughness.  The  roughness  measures  we  will  consider  in  most  detail  are  defined  by 

< 1  - 1 )  J, (u)  ■  /  (u  2  ♦  u  ‘)dx,dx_ 

ft  1  2  1 

or 


J2(u) 


.,2  2 

*  2u  ♦  u  )dx,dx. 
X.X-  x_x_,  1  2 

1  *.  £  £ 


We  will  also  briefly  consider  the  more  general  measures 


m  /  m  \  2 
T  (u)  -  f  l  (m)[ — j — J  dx  < 

’  o  iio  1  taX  '  1 


m  -  1,2,3, 


First  we  consider  the  problem  of  estimating  u.  With  the  roughness  measures 
(1.1)  and  (1.2)  our  estimate  u(Xj,x  )  of  ulx^.x,)  will  be  the  solution  to  one 
of  the  following: 

Problems  I-1/I-2:  Find  u  c  X  (an  appropriate  space  of  functions  on  ft) 
to  minimize  vTj(u)/J^(u)  subject  to  the  volume-matching  constraints: 

(1.4)  /  f  u (x , y) dxdy  ■  s  ,  i  *  l,2,«»»fN  , 

Qi 

N 

where  U  ft  ■  ft. 
i-1  1 

We  obtain  a  characterization  of  the  solution  to  a  general  problem  of  which  Problems 
1-1  and  1-2  are  special  cases.  This  is 
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p»  lem  1-4:  Let  4  be  a  smooth  bounded  subset  of  R  ,  Euclidean  d-space. 
Find  u  t  nm(Q)  to  minimize  J(u)  -  4(u,u),  where 


4(u,v)  -  y  /a  „Dau  D^v  dx 

|a| , | 8|«m  ft  U 


subject  to 


/  $  (x)u(x)dx  -  s  ,  i  ■  1 ,2,  •  •  •  ,N 

n  1  1 


Here  H  (ft)  is  the  Sobolev  space  of  functions  with  mixed  partial  derivatives  up 


to  order  m  in  L-(ft),  x  -  (x, ,x  , • • • ,x  .)  ,  a  ■  (a, .)  ,  B  ■  (B,, •••,&.), 

2  12d  Id  Id 

d  d  d  a  3a 

ja|  -  V  a,  7  i  -  y  B.  -  m,  pau  *  - - - ;  a  „  are  functions  of  x 

•i  1  •  .  1  i,  i  a,  a. 

i-l  i-l  i-1  J*  1. --Sx^ 

satisfying  certain  conditions  specified  in  Section  2,  and  the  (4^)  are  linearly 

independent  functions  in 


The  characterization  of  the  solution  to  Problem  1-4  is  qiven  in  Section  2. 
Certain  further  details  are  carried  out  in  Section  3  for  the  special  cases  of 
Problems  I-l  and  1-2.  A  simple  example  of  Problem  I-l  with  concentric  circles 
as  subdomains  is  worked  out  explicitly  in  Section  4. 

Numerical  algorithms  for  computing  the  solutions  to  Problems  I-l  and  1-2  will 
appear  in  a  separate  paper. 

The  solutions  to  problems  I-l,  1-2  and  1-4  are  not  required  to  be  non-neqative 
althouqh  it  is  known,  of  .course,  that  ulx^.x^)  and  glx^.Xj)  are  non-negative. 

In  this  paper  we  sidestep  the  philosophical,  theoretical  and  computational  problems 
of  imposing  non-negativity  on  the  solution,  and  hope  to  address  this  problem 
separately.  The  results  of  Lions  and  Stampacchia  (12)  will  be  relevant. 

We  know  of  a  very  short  literature  specifically  on  the  volume  matchinq  problem 
(Although  it  is  of  course  only  a  special  case  of  the  well  studied  problem  of  esti¬ 
mating  a  function  given  the  values  of  some  linear  functionals,  see  Golomb  and 


Weinberger  (91,  Kimeldorf  and  Wahba  (111.)  Boneva,  Kendall  and  Stefanov  (2)  dis¬ 
cuss  a  special  case  in  one  dimension.  Schoenberq  and  de  Boor  (16)  discuss  a 
volume  matching  problem  in  two  dimensions  where  the  roughness  measure  has  a  tensor 
product  structure  and  il  is  a  rectanqle  with  the  ft. 's  a  rectangular  subdivision. 

Our  interest  in  this  problem  was  sparked  by  a  paper  of  Tobler  |18].  He  proposed 

2  2 

to  solve  the  volume  matching  problem  by  minimizing  J.  (u)  *  f  f  ( u  ‘  +  u  ‘)dx  dx 

ft  X1  X2 

subject  to  volume  matching  conditions,  positivity  constraints,  and  certain  boundary 
conditions,  and  suggested  a  numerical  algorithm  for  doing  this.  Some  of  the  results 
here  are  alluded  to  in  our  comments  to  his  paper  (Dyn,  Wahba  and  Wonq  |8]). 

Our  results  show  that  the  solution  to  problem  I -A  and  the  special  cases  1-1 
and  1-2  satisfies  a  certain  eliptic  boundary  value  problem  with  Neumann  boundary 
conditions.  Numerical  implementation  of  these  boundary  value  problems  can  be 

avoided  if  one  is  willing  to  modify  the  roughness  criteria.  Let  X  be  a  suitable 

2 

space  of  functions  on  R  (to  be  defined),  and  define  J  on  X  by 


J  (u) 
m 


!  I 


i  '"’l 

i«0  \ 


,m 

3  u 


ax^Xj  ) 


dxidx; 


We  consider 

Problem  I-m:  Find  u  c  X  to  minimize  J  (u)  subject  to 

w 


/  u  dx  *  s.,  i  *  1,2,***,N 
°i 

If  u  is  the  solution  to  this  problem,  we  will  have 

J  (u)  >  J  (u)  >  J  (u)  , 

in  —  m  —  ro 

with  inequalities  holding  in  general.  This  approach  of  using  J  (u)  as  a  rough- 

w 

ness  criteria  has  been  extensively  used  for  estimating  surfaces  given  evaluation 
data  by  Duchon  |6] ,  (7],  Meinquet  (131,  Paihua  and  Utreras  (15)  and  Wahba  (19]. 
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we  derive  in  Section  7  an  explicit  ex 


Using  these  available  results 


pression  for  the  solution  of  Problem  I-m,  and  a  readily  computable  approximate 


solution.  The  results  generalize  easily  to  d  dimensions 


We  now  proceed  to  the  problem  of  estimating  g.  Since  the  data  Z.  are 


only  estimates  of  the  u.  we  only  want  g  to  satisfy  volume-matching  conditions 


approximately.  As  in  the  case  of  smoothing  splines  (see  { 5 ]  and  references  therein) 


Problem  I I-m:  Find  g  <  X  to  minimize 


(Z.  -  /  g(x,y)dxdy)  ♦  X  J  (g) 

1  *  nl 


with  J  (u)  defined  by  (1.3).  Here  the  {w  }  should  be  equal  to  1/variance  Z 


The  parameter  X  represents  a  tradeoff  between  the  roughness  of.  g  and  the 


which  is,  of  course 


infidelity  of  g  to  the  data.  The  variance  of  Z.  is  u 


unknown.  In  practice,  the  w  would  have  to  be  chosen  iteratively.  One  could 


The  resulting  esti 


since 


mate  of  g  is  then  used  to  get 


for  a  second  estimate,  etc 


In  Section  5  we  characterize  the  solution  to  problem  II  for  J  given  by 


In  Section  6  we  indicate  how  X  may  be  chosen 


to  approximately  minimize  the  predictive  mean  square  error.  In  Section  7  we  give 


an  explicit  representation  for  the  solution  to  Problem  Il-m  with  J  replaced  by 


(Problem  Il-m) .  More  significantly,  we  give  explicit  formulae  for  approximate 


solutions  to  Problem  Il-m  which  are  suitable  for  numerical  computation.  In  this 


context  we  also  derive  formulae  for  computing  an  optimal  X  based  on  the  results 


Hopefully,  these  results  will  provide  the  first  step  towards  efficient 


methods  for  converting  aggregate  data  to  density  maps 


2.  Smooth  Surfaces  On  Bounded  Domains  Matching  Inteqral  Data. 

Consider  a  bounded  domain  of  with  T  its  boundary,  and  a  bilinear 


'Uu.u)  ■  )  |  a  lx)DUuD‘v  ,  a  (x)  *  L  (13) 

|«|.|B|-m  fl  “8 


where  x  »  (x^, • • • .x^) ,  a  »  (a ^ , • • • .a^) ,  |a|  *  j  a^,  a^-non-negat ive  inteqer, 

_  i-1 

,  1  ,  d 

ci  3  3 

D  "  -  •••  -  (and  similar  notations  for  6).  With  this  definition  /4(u,v) 

“l  .  “d 

3x,  3x. 

l  a 

is  continuous  on  where  H171  { Q)  is  the  Hilbert  space 

-  {u  !  D®u  t  L2(fl).  'o|  in),  Hull2  -  [  II  Dku  II2 

Hm(Q)  |k|-m  L2  (13) 

By  assuming  that 


1  »  <x) y  y  >  C  l 

_ ao  «»  8  0 


' a | ,  fi-m 


i  ci  |  «m 


for  all  y  »  (yi»,**,y<),  k  •  ■  ■  a  i  »  m  ■ ,  M(u.u))1  "  is  a  seminorm  on  Hm(fi) 

with  a  null  space  0  -  the  space  of  all  polynomials  of  total  degree  less  than  m, 

which  is  of  dimension  M  5  (P+<!  S. 

d 

In  this  section  wo  prove  the  existence  and  uniaueness  of  the  solution  to 
Problem  1-4: 

For  given  s,,’**,sf.,  find  u  t  Hn(fi)  minimizing  4(u,u)  among  all  functions 
m 

in  H  (13)  satisfying  the  integral  data 

(2.3)  /  u+ .  »  a  ,  i  *  1,***,N 

13  1 

2 

where  ♦,•••, ♦  are  N  linearly  independent  functions  in  L  ((3). 

In  particular  we  characterize  the  solution  of  Problem  1-4  as  a  solution  of  a 
certain  boundary  value  problem. 
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We  ve  two  lemmas. 


,  ::ima  2.1;  In  the  subspace  H ^  of  H  (fl)  given  by 


H  ■  {u  1  u  (  Hm(fl)  ,  /  d"u  »  0,  1  ai j  <  m) 

0  n 


/I(u.u>  is  an  equivalent  norm  to  llu 


hV:) 


Proof :  By  (2.1)  there  exists  >  0  such  that 


A  (u,u)  <  C.llull  ,  u  t  Hm(ft) 

"  1  Hm(n) 


Iterating  the  Poincare  inequality  (14] 


/  u2  <  C{  I  /  <Dau)2  ♦  (/  u] 2 } #  u  €  H‘<fl) 

a  ~  ]  a  |  *  i  n  n 


we  obtain  for  any  0  ^  k  <  m 


(2.7)  l  j  (Dau)2  <C.{  l  j  (Dau)2  *  l  if  Dau]2),  u  <  Hm(n) 


|a'-X  a 


!oj*m  fi  k^Joj<m  0 


Thus  by  (2.4)  and  (2.2) 

2  -  2  C1 

(2.8)  Hul!  <  C,  T  /  (Dau)  <  — -  A  (u,u)  ,  u  c  H  (D) 

"  3  i  a|*tn  0  "  C0 

bet  Q  •  span'! q.  ,  •  •  *  ,q  }.  We  assume  that  N  >  M  and  that  the  N  linear 
1  M 

functionals  in  (2.3)  are  linearly  independent  over  O.  Without  loss  of  generality 


we  can  assume  that  the  matrix 


qi*N-M+i1i,j-l 


is  of  rank  M.  Therefore  there  exists  a  basis  {q1# of  Q  with  the  pro¬ 


perty 
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(2.10) 


q  A  -6..,  i,j-l,***,M 

ivN-M+j  13  J 


Lemma  2.2:  In  the  subspace  of  H  given  by 


(2.11)  H  cm  -  {u  |  u  €  H®(R),  /  u*  -  0,  i  -  N-M-H,  •••  ,N) 

1  R 


•  4(u,u)  is  an  equivalent  norm  to  Hull 


Hm(fl) 


Proof :  For  any  u  t  HjlR)  there  exists  q  f  0  such  that  u()  »  u-q  <  H0(R) , 

and  therefore 

M 


u  •  U0  +  l  uoVM+i 

l-l  (1 


Since  for  any  $  <  L.}(R) 


(2.12) 


j/  u*|  <  Hull  ,  u  <  H  (R)  , 

L  (fl)  Hm(fl) 


fl 


we  get  in  view  of  Lemma  2.1 


luH  <  llu  II  ♦  T  II q.  II  IU  M  .  II  ,  Hu.H 

»"«!>  '  °  h"(Q)  1=1  1  N-M*‘  l2(OT  0  «■(!» 


<  C. ’A (u  ,u  )  -  C.*a (u,u) 
—  4  0  0  4 


This  together  with  (2.5)  completes  the  proof  of  the  lemma. 


(2.13) 


Let  u  (  H  (R)  satisfy  (2.3).  Then 

M 

“  "  U  '  l  SN-M.i  *i  C  ”l 


i-1 


(2.14) 


n 

(/*)  ■  •)  -  j,  Vh.i  /  ' Vj  *  V  1 


«  1 , • • • ,N-M  , 


and  A  (u,u)  -  >t(u,u).  Therefore  Problem  l-A  is  equivalent  to  Problem  (I-^)': 


Find  u  (  Kj  minimising  /Wu,u)  among  all  functions  of  Hj  satisfying 


(2.14),  or  equivalently  satisfying 
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(2.15) 


j  -  1  ,  *  •  • , N-M 


i4i  '  V 


with 

(2.16) 


»  *  -  y  a  4  ,(a..l  arbitrary 

j  ’j  ii1  i3VM+r  ij 


In  particular  it  is  possible  by  assumption  (2.9)  to  choose  )  such  that 


(2.17) 


f  qi  m  0,  j  ■  1 ,  •  •  • ,  N-M ,  q  f  0 

n  5 


Ry  Lemma  2.2  the  linear  functionals 


(2.18) 


L.(u)  •  J  u$..,  j  ■  !,•••, N-M 


1/2 


are  bounded  in  H ^  with  respect  to  the  norm  |4(u.u) J .  Invoking  the  Riesz 
representation  theorem  we  conclude  the  existence  of  ,  j  *  1»***,N-M 

satisfying 


(2.19) 

and  due  to  (2.17) 

(2.20) 


/4(u,C  )  "  /  u4  .  all  u  <  Hj 
0  3 


/4(q,C.)  -  /  q$.  "  0,  all  q  <  Q 

3  n  3 


Since  are  linearly  independent  so  are  i  •  *  *  * ' ^n-M'  and  t*le  so*ut*on  to 

Problem  (I-/1)'  is  known  to  be  the  unique  function  in  the  span  of 
satisfying  (2.16),  (see  [9)).  The  solution  to  Problem  I -A  is  related  to  this 
solution  according  to  (2.13).  The  following  theorem  summarizes  the  above  findings: 

Theorem  2.1:  There  exists  a  unique  solution  to  Problem  I -A.  The  solution  is  of 
the  form: 


(2.21) 


N-M  M 

”  ’  CiEi  4  A  %-M*i*i 


i-1 


i-1 


■9- 


f 


•f*frr 

where  are  the  unique  functions  in  H,  determined  by  (2.19),  and 

•••#«?  are  the  solution  of  the  non-sinqular  linear  system 


N-M 


W-"’  Jx  •*)•*)■  j,  *-■*,  / 


Vj*  j  "  ' 


An  immediate  consequence  of  Theorem  2.1,  (2.19)  and  (2.20)  is 
Corol lary  2.1:  The  solution  u  of  Problem  1-4  is  uniquely  determined  by  the 
variational  characterisation 

N 


(2.23) 


4(u,v)  -  /  (  1  y j».)v,  v  <  H®(R) 
n  i-i 


and  the  matching  conditions 


(2.24) 


i  -  1  ,  •  •  •  ,N 


In  (2.23)  y  are  constants,  which  in  particular  satisfy 

1  « 


(2.2S1 


J  <  l  Vi)q  *  °«  9*0 

Q  i-1 


In  case  fl  is  a  smooth  domain  the  solution  u  of  Problem  1-4  can  bo 
further  characterized  in  terms  of  a  boundary  value  problem.  Since  each 
1  <  i  <  N-M,  satisfies  (2.19)  and  (2.20),  namely 

4(u,t  )  *  /  u4  for  all  u  <  HB(R)  , 

R 


we  conclude  from  Corollary  2-2  on  pages  219-220  of  Aubin’s  book  (l]  that  C. 
is  the  unique  solution  in  to  the  boundary  value  problem: 


(2.26) 


AC 1  ■  in  R  , 


(2.27) 


5  Ct  »  0  for  n  <  )  <  2m- 1  on  T  . 


In  (2.26)  A  is  the  differential  operator  of  order  2m  given  by 
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(2.28) 


and  in  (2.27)  6 

mapping 


(2.2 9) 


r  I  fil  8  a 

Au  -  £  (-1)  1  'd  (a  ^(xlD  u) 

|a| . [B|-» 

( c  .,•••,*  )  is  i  differential  operator  of  order  m 
2m- 1  m 

HB((1.A)  ■  (u  u  <  HW(Q)  ,  Au  (  t.'  (.0  ) 


1 

m  m-)-  - 

into  n  H  (D,  such  that  the  generalized  Green  Formula  holds: 

j»2m-l 

<n;l  /3j  \ 

(2.30)  <4(u.v)  -  /  ( Au) v  ♦  /  /  (6  u)l — rvj 

Q  j-0  Q  2B"3'1  W  / 

(—  is  the  operator  of  normal  derivative  to  the  boundary  "). 


The  characterization  (2.26),  (2.27)  of  (^  •  toqether  with  Theorem  2.1 

yields 

Theorem  2.2:  The  •  1  it  n  to  Problem  Z-A,  for  a  smooth  domain  0,  is  uniquely 
determined  as  the  solution  to  the  boundary  value  problem 

N 

(2.31)  Au  •  l  Vi*i  in  fi 


(2.32)  i  u  •  0  m  *  J  <  2m_1  on  r 

which  satisfies  the  matching  conditions  (2.24).  Tn  (2.31)  are 

N  constants  satisfying  (2.25)  . 


-11- 


3.  [.aplacian  Histosplines  -  The  Volume-Matching  Surfaces. 

In  this  section  we  specialize  to  the  concrete  problem  of  findinq  a  smooth 

2 

surface  u  -  u(x},x,)  having  prescribed  volumes  over  specified  subdomains  in  R 
We  characterize  the  volume  matching  surface  as  a  function  with  the  even  order 


differential  form  A 


*  IM*  I  AY 

•  Ikj  *  NJ 


of  constant  value  in  each  of  the  sub- 


domains.  These  surfaces  are  therefore  strikingly  analogous  to  even  degree  one- 
dimensional  splines,  regarded  as  functions  with  a  certain  even  order  derivative 
of  constant  value  in  each  subinterval.  Following  a  suggestion  of  Professor  Iso 
Schoenberg  we  term  these  surfaces  "Laplacian  Histosplines”  in  analogy  to  the  uni¬ 
variate  Histosplines  of  Boneva,  Kendall  and  Stefanov  (21,  which  are  the  even  degree 
univariate  splines  solving  the  "area  matching”  problem. 

We  consider  in  details  the  following  two  problems:  Let  £1  be  a  smooth 

2  N 

bounded  domain  in  R  subdivided  into  N  disjoint  domains  0,  *•••,0  ,  £1  -  <->  fl.  . 

1  N  .  ,  i 

i-i 

Problem  I- 1 :  Find  u  <  H^(£l)  minimizing 


/  (u  '  ♦  u  “)dxdx. 

,  x.  X,  l  . 


Q  1 


among  all  functions  in  H  (ft)  satisfying 


/  u  -  s  ,  1  -  1 • • • • »N  . 

Qi 


Problem  1-2:  Find  u  <  H  (Q)  minimizing 


/  (u‘  ♦  2u^  ♦  u‘  )dx,dx 

x,x,  x,x„  x^x„  1 


n  ii 


1  2  2  2 


amonq  all  functions  in  H  (£1)  satisfying  (3.2). 

From  a  practical  point  of  view  these  two  problems  are  the  most  interesting, 
since  computation  of  solutions  of  similar  problems  with  higher  order  forms  (2.1) 
becomes  too  complicated,  with  the  increased  complexity  of  the  operators  A  and  6 


in  Theorem  2.2. 
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Usinq  Theorem  2.2  for  the  special  sett inq  of  Problem  1-1  toqether  with  the 
classical  Green  Formula  [31: 


fuv  ♦  u  v  -  /  (-Au)v  ♦  /  v 
4  x  x  x  x  *  '  3n 

Q  1  1  2  2  Q  r 


we  obtain: 

Theorem  3.1:  The  solution  to  Problem  1-1  is  uniquely  determined  by  the  following 
conditions : 


i-l  i  i 


1  in  fl . 


0  elsewhere 


I  yj  i-o 

i-i  fij 


0  on  r 


/  u  •  s  ,  i  -  1  #  •  •  •  «N 

Qi 


To  get  a  similar  result  for  Problem  1-2,  we  first  derive  a  more  general  Green 
Formula  for  the  bilinear  form  corresponding  to  the  seminorm  (3.3).  By  a  repeated 
use  of  (3.4)  we  get 

0.5)  /  u  v  ♦2u  v  +u  v  -  /  (A2u)v  -  /  {-£-  Au)  v  ♦  /  V  •  Vv 

a  *1  X1  *i*2  X1X2  X2X2  X2X2  n  r  *n  r  8n 

since  on  T  ?u*Vv  -  r—  r—  ♦  •1^  r— ■  ,  where  —  is  the  tanqential  derivative  to 

Jn  dn  3t  3t  3t 

f,  the  last  term  in  (3.5)  becomes 


/HS 


riH.w  .  r  lly.  tL  4  r  ji L  .  f  ii i 

r  3n2  Jr  Jr  ^ 


32u  3v  m  r  3^u  3v  r  v 

3x3n  3t  Ij.  3n2  3n  Jt23n 


Comparing  (3.5)  and  (3.6)  with  the  generalized  Green  Formula  (2.30),  we  conclude 
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that  for  the  seminorm  (3.3),  A  and  6  of  Theorem  2.2  are: 


(3.7) 


»-«•-.  s  -  ,vv.  i 


Thus  by  Theorem  2.2: 

Theorem  3.2.  The  solution  u  to  Problem  1-2  is  uniquely  determined  by  the 
following  conditions: 


.2- 


“  ■  .1,  yi  *0,  ln  0 


i-1 


N 


l  /  q  -  0.  a  -  1.x., X 
i-1  (1 


-  0,  (-“  ♦  A)  —  -  0  on  T 

3n2  3r2  3n 


Q 


u  «  s  ^ ,  i  N 


i 


Remark:  It  can  be  shown  by  Theorem  2.2  and  repeated  applications  of  the  classics 

Green  formula  that  for  the  higher  order  roughness  criteria 

2 

m  /  _m  i 

r  r  m  I  3  u 

(3.8) 


Vu)-/  jOy-rirr)  dKidv 

n  i-1  \3x13x2  / 


the  solution  to  the  volume  matching  problem  satisfies 


(3.9) 


(-1)  A  u  -  l  y.Xn  In  0 
i-1  i 


with  appropriate  boundary  conditions 


(3.10) 


fl^u-0  on  T  m  <  j  <  2m- 1 


14- 


The  total  number  of  linear  equations  (4.3)  -  (4.7)  is  3n,  as  is  the  total 
number  of  unknown  coefficients.  If  R(J  >  0  there  is  an  additional  boundary  con¬ 
dition  to  be  satisfied 


(4.8) 


du 

dr 


r-R 


Vo  ci 

—  if  R0>0  . 


Claim:  If  Rq  >  0,  (4.8)  is  linearly  dependent  on  equations  (4.3)  -  (4.5). 

R.  ■  0  then  (4.3)  -  (4.5)  imply  c^  •  0. 


If 


Proof :  Summinq  (4.5)  for  i  -  1,**-,N-1  we  get 
N-l  R2  ,  N 

y  (y  ~y  )_i  -  I  y 

t  i  +  l  2  2  .  *• 

i»l  l-l 


ii  a  r»  .  ,  _  _  _  _  _  _ 

r  .  l  1  r  2  2,1  _2  1  2 

C1  "  "N  "  <>i"Yi  +  l)  2  “  2  Yi  Ri~Ri-l  *  2  Vo  "  2  >NRN 


which  in  view  of  (4.4)  and  (4.3)  can  be  written  as 


Ci  -  CN  '  2(Y1B0“YNRN>  "  2  Vo  ‘  CN 


Therefore  c^  -  j  y^R*,  proving  the  claim. 

By  integrating  ru(r) ,  one  can  transform  this  volume  matchinq  problem  into 
an  interpolat ion  problem.  (Similar  equivalence  exists  between  area-matching 
splines  and  interpolat inq  splines  in  the  one-dimensional  case  (16)).  Thus  defining 


(4.9) 


U(r)  -  /  ou(o)do  ,  u(r)  -  —  U'(r) 

Ro 


we  have  to  construct  an  "interpolating  spline"  of  the  form: 

2  4  2 

(4.10)  U(r)  -  ♦  B^r  ♦  c^r  4  D^r  log  r,  R._^  <  r  <  Ri ,  i  -  1,***,N 

satisfying 


(4.11)  U(r)  *  C  (R  .R  )  , 

0  N 


i  k 

U(R.)  -  f  s.  ,  i  -  1»  *  * • *N  • 

1  2w  J-l  1 


16- 


— 


2  4  2 

It  is  easy  to  check  that  the  functions  I,  r  ,  r  ,  r  log  r  constitute  an 
Kxtended-Chebychev-System  on  any  interval  of  the  form  (0#R^).  Thus  U(r) ,  con¬ 
sidered  as  a  function  of  r,  is  a  Chebychev-spline.  (For  the  notion  and  construct 
ion  of  Chebychev  splines  see  e.g.  (10)  Chapter  10) . 


5.  Laplacian  Histosplines  for  Inexact  Data. 

In  this  section  we  consider  the  problem  of  finding  a  smooth  function  g 
given  inexact  volume  data.  Similar  analysis  can  be  done  in  the  more  general 
settinq  of  Section  2. 

Problem  II-m-.  For  a  given  set  of  data  zi'***'2N  find  g  e  Hm(fi)  mini 
mizing 

N 

(5.1)  l  w  [J  g  -  Z  )  ♦  XJ  (a) 

i-1  1  fi.  “ 

i 

where  J  (g)  is  defined  in  (3.8),  fi,  fi,  ,•••,(!„  are  as  *.n  Section  2  and 
w  IN 

X,  w  ,•••,*  are  fixed  positive  constants. 

1  N 

In  the  notation  of  Section  2  any  g  c  Hm(fi)  can  be  represented  as 
9  "  ql  4  g2  4  93  where  g^  f  Q,  t  span{ Cj . • • • . ^N_M)  and  q ^  satisfies 

(5.2)  J  g  -  0,  i  -  1,*«*,N 

ft  ^ 


By  (5.2)  g  c  H  is  orthogonal  to  £,*••,£  „  with  respect  to  the  inner- 

•J  1  1  N—M 

product  in  H,  corresponding  to  the  norm  » J  ( • ) .  Therefore  g,  does  not 
1  m  3 

affect  the  first  term  in  (5.1)  while 


V91  4  9 2  4  93) 


W 


4  W 


and  necessarily  the  solution  to  Problem  II-m  is  of  the  form 

N-M  M 


(5.3) 


9  *  9X  ♦  <?2 


l  ci^i  4  l  di9i  • 
i-1  i-1 


Since  for  the  volume  data,  in  Section  2  are  of  the  form 

1  N 


♦i  "  *n/  1  "  ' 


hence  by  (2.16),  (2.17)  and  (2.26) 


(l 


(5.4) 


ci  •  \  0 


in  Cl. 


in  n,  j  ^  i,  j  ■  1 , • • • , N-M 

j 

in  ft ,  j  »  N-M+l , • • • ,N 


(jiJ  5 


with  y  satisfyinq 


N 

(5.5)  V  >  ./  q{  ♦  /  q.  ■  0,  t  -  1,***,M,  i  » 

j-N-M+1  13 SI.  ci, 

3  l 


In  view  of  (5.4),  (5.5)  and  (2.27)  the  solution  ci  to  Problem  11-m,  qiven  by 
(5.3),  satisfies  the  boundary  value  problem: 

N 


(5.6) 


(-1)  A  q  -  l  y.xQ  in  Cl 
i»l  i 


(5.7) 


i  ,  q  *  0  on  T  m  <_  j  2m- 1 


with  >  ,••*,>  N  constants  restricted  by 
1  N 


(5.8) 


l  Yj/  q  .•  0,  l  •  l.—  .H 

i-i  n. 


In  (5.7)  the  boundary  operators  6  .*••,6.  ,  are  as  in  the  Remark  in 

m  2m- 1 

Section  3. 

The  followina  theorem  relates  the  values  of  the  constants  Y, .•••»>„  in 

1  N 

(5.6)  to  the  "smoothed  data”,  namely  to  the  values 

(5.9)  Z  *  /  q,  i  -  1  * • • • ,N  . 

ni 

Theorem  5.1.  The  solution  q  of  Problem  Il-m  satisfies  (5.6)  with 


(5.10) 


-  — (Zi-Zi),  i  »  1 , • • *  «N  . 
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Proof :  The  coefficients  in  (5.3)  satisfy  the  necessary  conditions  for  minimizing 

(5.1),  namely  the  vanishing  of  the  partial  derivatives  of  (5.1)  with  respect  to 

c  , •••,c  „  and  d  .•••,d„.  In  terns  of  the  bilinear  form  A  (•,*)  correspond- 

1  N-M  1  M  m 


r  '  m 


ing  to  J  (•),  these  conditions  become: 
m 


(5.11) 


l  'VV5 '/  S  *  *VWei> 

i-i  n. 


o,  j  -  i , • • • ,N-M  , 


(5.12) 


l  W  (Z  -Z  )f  q.  -  0,  j  -  1 , • • • *M 
i-N-M+i  n  J 


In  deriving  (5.11)  wo  recalled  that 


(5.13) 


/  C.  ■  0,  i  ■  N-M+ 1 , • • • , N ,  j  ■  1 »  * • • , N-M 

£1.  3 

i 


Let  K  bo  the  (N-M)  *  (N-M)  matrix  with  entries 


(5.14) 


!d  ■  wv  ■  L  cj  •  /  v 


-  1 , • • • ,N-M  , 


let  T  be  the  (N-M)  *  M  matrix  with  entries 


ru  ■  L 


,  i  -  1 ,  • • • ,N-M,  j  ■  1 ,  • • • ,M 


and  let 


W  -  diag{w1,...,wN_M),  W  -  diaq{wN_M+1,...,wN)  . 

c  -  <ci'*“'cn-mv'  2  ■  ‘V  —  'Vm**  ; "  ‘Wr* 


'V  ' 


1  ?  *  (ZN-M*1',”'ZN)' 


With  these  notations  (5.11)  and  (5.12)  become 


(5.15) 


KW(z-z)  -  iKc  -  0 


(5.16) 


z  -  z  *  -w"1T'*f(z-i) 
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Since  K  as  defined  in  (5.14)  is  symmetric  positive  definite,  (5.15)  implies 


C  *  —  W ( 2-2) 


(5.17) 


while  by  (5.3),  (5.4)  and  (5.6) 


(5.18)  ■  (-l)mAmg  ■  y.  in  fV  ,  i  -  !#•••, N-M 


Therefore  (5.10)  holds  for  i  «  !,•••, N-M,  and  (5.8)  becomes 


( S . 1 9)  ( v  #•••*>)'■  -T '  ( Y  » • • • • Y  )' 

1  '  N-M* i  N*  VY1#  'TN-m' 


-  T*c  ■  -  yT’W(z-z) 


Comparing  (5.19)  with  (5.16)  we  conclude  that  (5.10)  holds  for  i  -  N-M+1,***,N 
as  well. 

A  direct  consequence  of  Theorem  5.1,  the  representation  (5.3)  of  g  and 

(5.4)  is: 

Corollary  5.1:  The  solution  of  Problem  Il-m  is  of  the  form 


(5.20) 


*  "  T  l  wi(zi"iiKi  4  *  Wi 

i-1  i-1 


and  satisfies  the  integro-differential  equation 


(5.21) 


(-l)Vg  -  j  l  Xn  w  [Z  -  /  g] 
i-1  i  1  n 


with  boundary  conditions 


(5.22) 


S^g  -  0  m  <_  j  <  2m- 1  . 


Equations  (5.21),  (5.22)  indicate  an  alternative  direct  way  for  the  computa¬ 
tion  of  g,  avoiding  the  computation  of  the  functions  £,•••,£ 

1  N-M 

We  conclude  this  section  by  deriving  explicitly  the  relation  between  the 

vector  of  given  data  Z  -  (Z  ,•••,2  )’  and  the  vector  of  smoothed  data 

1  N 

z  -  (Z^....^)*. 
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From  (S.20)  we  get 


z  -  Y  KW(z-z)  ♦  7* 

and  after  substituting  for  z  from  (5.16) 

-  -1  *  1 

(5.23)  z  -  Tz  ♦  TW  T’W(z-z)  ♦  r  KW(z-z) 

With 

B  -  (I  +  j  KW  ♦  Tw’Vw)"1 


(5.23)  and  (5.16)  become 


(5.24) 


Z  -  z  ■  B(Z-Tt),  z  -  z  -  -W  ^T'WBfZ-Tz) 


Combining  the  last  two  expressions  we  conclude  that 


(5.25) 
with 

(5.26) 


Z  -  A  (  X)  Z 


I  -  A ( X) 


-BT 


\-V~ ^T' WB  W-1T'WBT / 
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6.  The  problem  of  choosing  X. 

We  give  a  procedure  for  choosing  X  in  Problem  II.  In  this  section  we  suppose 
(inaccurately!)  that  the  { w ^ )  in  the  definition  of  Problem  II  are  given  positive 
constants.  In  the  problem  presented  in  the  introduction  we  want  w^  -  1/variance  Z^ 
*  1/u^ .  Since  the  r  ^  are  being  estimated,  the  w^  can  be  chosen  iteratively 
by  one  of  several  obvious  ad  hoc  procedures.  In  what  follows,  the  w^  are  assumed 
fixed  and  given.  It  is  likely  that  w^  1  will  give  reasonable  answers  in  most 

cases  when  the  i.  are  all  of  the  same  order  of  magnitude. 

A  good  criteria  for  choosing  X  is  the  minimization  of  R(X)  defined  by 

(6.1)  R(X)  -  E  l  6  (p  -  /  q  )2 

i-1  11  Jl. 

i 

where  E  is  expected  value  and  g^  is  the  solution  to  Problem  II  and  the  ©^ 
are  given  positive  weights.  Since  the  u ^  are  not  known,  we  cannot  minimize 
R(X).  However,  an  unbiased  estimate  R(X)  of  R(X)  is  available  by  generalizing 
an  observation  in  Craven  and  Wahba  (5).  Let  A(X)  be  the  N  *  N  matrix  satis¬ 
fying 


Such  a  matrix  is  given  explicitly  in  (5.23),  (5.24). 

Then  (6.1)  becomes 

R(X)  -  E  Hd*/2(u  -  A  (  X)  Z)  II  2 

where  D  *  diag{0.  ,*•*,©) ,  and  u  *  (r  ,  • • •  ,u  ) '  .  Defining  e  -  (  ?.,•••,€)'  by 


z 


u  *  t 


we  have 


[(I  -  A  ( A)  )  Vi  -  A  (  A ) «  H 


(U  -  A  (  A )  Z )  II 


(I  -  A<  A)  vjII  ♦  Tr  D  A  (  A)  £  A  *  (  A ) 


where  E  ■  diag(var  Z  .  var  Z 


Z  ! .  We  claim  that  an  unbiased  estimate  R ( A )  of  R(X) 


Let  Z  •  diagi  Z 


is  given  by 


(I  -  A  ( A) )  £  ( I  -  A  (  X )' )  D, 


(I  -  A  ( A ) )  Zll 


A  {  X)  E  A  (  X D 


In  fact  (6.2)  Simplifies  to 


0.Z,  -  2  Trace  D.KI  -  A( X) ) 

11  O 


(I  -  A ( X) ) Z B  ♦ 


To  assert  our  claim  observe  that 


(I  -  A ( A ) wH 


(I  -  A( X ) ) Z 


(I  -  A( X) ) E ( I  -  A(X) ' ) D 


Substituting  (6.4)  into  (6.2)  and  using  (6.5)  we  obtain 


Thus  it  is  reasonable  to  choose  X  by  minimizing  R(X) 


7.  Laplacian  Histosplines  for  a  Modified  Smoothness  Criteria. 

Problems  in  coding  a  numerical  algorithm  for  computing  u  and  g  related 
to  solving  the  Neumann  boundary  value  problem  in  irregular  domain  can  be  avoided 
by  modifying  the  smoothing  criteria  somewhat. 

Whether  or  not  this  modified  smoothing  criteria  gives  results  equally  pleasing 
as  the  smoothing  criteria  previously  used,  and  whether  the  computing  time  required 
is  comparable  or  not  remain  to  be  seen.  However,  the  coding  of  an  algorithm  for 
the  modified  criteria  appears  to  be  relatively  straightforward,  and  is  similar  to 
already  oxistinq  codes  for  the  case  of  point  evaluation  data  (13),  [151,  (19). 

The  results  below  are  modest  generalizations  of  results  given  by  Duchon  [6], 
[7],  and  later  discussed  by  Meinguet  [13)  and  Wahba  [19]. 

We  let  d  ■  2,  however  the  generalization  to  arbitrary  d  dimensions  is 
immediate  from  the  known  results  whenever  2m-d  >  0.  Let  X  be  a  suitable^ 
space  of  functions  on  R*  for  which 


J  (u) 
m 


l  'Kife)’ 


is  well  defined  and  finite. 

We  modify  problems  1-m  and  Il-m  to  the  following: 

Problem  I-m:  Find  u  t  X  to  minimize  J  (u)  subject  to 

in 

/  u(x.  ,JC  )dx  dx-  -  S  ,  i  -  1 , 2 , •  •  •  ,N 
Qi 

Problem  Il-m:  Find  g  «  X  to  minimize 


l  *^[2*  -  /  gfx^x^dx^Xj]  +  X  3B(g) 
i-1  Q. 


X  is  the  Beppo  Levi  space  of  all  the  Schwartz  distributions  for  which  all  the 
partial  derivatives  in  the  distributional  sense  of  total  order  m  are  square  in- 
teg  r  able  in  R2  (13)  . 
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whore  a.  *  |$V|  being  the  area  of  E.  and  n^  the  number  of  mesh 

points  in  Q. .  We  now  consider 

Problem  I-m-{t  i:  Find  u  «  X  to  minimize  J  (u)  subject  to 

"  m 


a  l  u  ( t . )  «  s  ,  i  -  1 , 2 ,  •  •  • ,  N  . 

t,<n.  1 

t  l 


Problem  II-m-{t  1 :  Find  g  >  X  to  minimize 


l  w  (Z  -  at  [  9(t.) 1'  ♦  X  J  (g) 

i-1  t,«Q. 

l  i 


Theorem  7.1:  Suppose  the  N  *  M  matrix  T  with 


(7.2) 


jv 


V°j 


W 


is  of  rank  M.  Then  the  solutions  to  problems  I-m-it,}  and  are 

unique  and  have  representations: 


(7.3) 


where 


and 


<qv(*)> 


M 

1 


span  the  space  of  polynomials  of  total  degree  less  than  m. 


T1 


coef f icients 


c 


(c 


r 


,cN> ’  and  d  - 


satisfy  the  followinq  equations: 
Problem  I-m-ft^l 

(7.4) 

(7.5) 


Kc  ♦  Td  ■  s 
T'c  ■  0 


where  K  is  the  N  *  N  matrix  with  ij  th  entry 


K.  .  -  a  a.  I  E  (t.  ,t.) 

lJ  1  1  i  i.  ■  k  1 
k  l 

t.€fl  . 

»  3 


i,j  ■  1 # •  •  •  *N 


and  s  “  (#_#•• • ,s) '  . 

1  N 

Problem  Il-m- ( t^  } 

(7.6)  (K  ♦  XW~ 1 ) c  ♦  Td  -  z 

(7.7)  T’c  -  0 

where  W  •  diaq (w  , • • • ,w  },  and  Z-  (Z. , • • • ,Z„) ' . 

IN  IN 

Proof:  The  special  case  n4  »  a  »  w.  «  1,  i  ■  1,2,***,N  is  just  the  problem  of 
interpolating  or  smoothing  evaluation  data,  and  in  this  case  the  result  has  been 
given  explicitly  in  (6],  (7J ,  (131,  (191.  The  extension  to  the  case  of  general 
n^,  a^  and  w^  is  straightforward  from  these  results  and  is  omitted. 

Observe  that  the  solution  to  problem  I-m-{tj}  can  be  obtained  by  solving 
(7.6)  and  (7.7)  for  the  solution  of  Problem  II-m-(tt),  with  X  -  0  and  Z  replaced 
by  s.  We  now  put  equations  (7.6)  and  (7.7)  in  a  form  suitable  for  the  computation 
of  c,  d  and  R(X).  Let  R  be  any  N  «  (N-M)  matrix  satisfying  R'T*0  . 

Since  T'c  *  0,  there  exists  a  unique  N-M  vector  b,  say,  with  c  ■  Rb.  Left 
multiplying  (7.6)  by  R*  and  substituting  c  *  Rb  gives 
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We  next,  assert  that  R'KR  is  strictly  positive  definite.  To  prove  this  we  use 


the  following  result  16] 


t  do  not  fall  on  a  straight  line.  Let  f  ■  (f 


any  non  2ero  vector  satisfying 


f  f,E  (t,-t.)  >0.  We  need  to  show  that  if  r  »  (r 


then 


satisfies  T'r  “  0,  then  r'Kr  >  0.  Let  F  be  the  n  *  N  matrix  with  jk 


and  0  otherwise,  let  E  be  the  n  *  n  matrix  with 


entry  a, 


entry  q  (t.) 


entry  E  (t  -t  ) 


and  let  T  be  the  n  *  M  matrix  with 


Then  K  -  F'EF  and  T  -  F'T.  Suppose  T'r  ■  0.  Then,  if  f  ■  F’r,  we  have 


In  case  X  *•  0  or  X  is  a  given  positive  constant,  b  is  obtained  from 


c  -  Rb  and  d  is  obtained  from  (7.6)  as  the  solution  of  the  system: 


pute  R(X)  we  first  obtain  an  expression  for  A(X).  The  appropriate  definition 


Using  the  fact  that 


one  obtains  from  (7.3) 

(7.10) 


*1**1 


yv 


*u 


A ( X)  Z  '  Kc  +  Td  . 


Combining  (7.6)  and  (7.10)  we  qet 

(I  -  A ( X ) ) Z  «  (K  ♦  XW_1)c  ♦  Td  -  (Kc  ♦  Td)  -  XW_1c  . 


Since  by  (7.8)  and  the  definition  of  b 

c  «  Rb  -  R(R'(K  ♦  XW_1)R)‘1R'Z 

we  finally  obtain 

(7.11)  I  -  A (X)  -  Xw"lR(R’(K  ♦  XW_1) R) ”^R*  . 

R  can  always  be  chosen  so  that  R'W  *R  •  I  ,  giving 

N*n 

I  -  A ( X)  -  XW_1R(B  ♦  X I) “ XR* 

where  B  »  R*KR  is  a  symmetric  positive  definite  matrix.  Now  let  UDBU' 
the  eigenvalue  decomposition  of  B  with  DB  *  diagf b^ # • • • #b^_^} >  then 

(7.12)  I  -  A( X)  -  XW_1RU(D  ♦  Xl)“1U’R'  . 

*3 

Recalling  the  expression  (6.3)  for  R(X): 

N 

(7.13)  R(  X)  -  II  Dq  z  (I-A(X) )  zai  ^  ♦  l  e.Zi  -  2  trace{DeI(I-A(X)) ) 
and  substituting  (7.12)  we  obtain 
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be 


(7.14) 


R(  X) 


*S 


N-M 


l 


l  h 

j*i 


i) 


b .  ♦  X 
1 


N  N-M 

y  o.z.  -  2x  y 

1=1  ,  1=1 


b.+X 

l 


where  v  *  (v  ,  •••,v  ‘  “  U'R'Z  , 

I  N“M 

0  6 

H  -  (h.  .}  -  U '  R '  W  1D.w”1RU  -  U'R*  d iaq {— i-,  •  •  • , RU 

II  o  it 

w  w 

1  N 

-1  Z161  Zn9N 

L  -  U.  .}  »  U’R'ED.W  Kl)  «  U'R*  diaq{  — — ,•••, - -}RU  . 

13  e  Wj  «n 

In  the  special  case  D.  *  W,  the  matrix  H  is  I  since  R'W  *R  =  I 

V 

(7.14)  simplifies  to 


N-M 


(7,15) 


R(  X )  -  X l 


i=l  (b.+X) 


7  w.z.  -  2X  y 

L ,  i  1  L 


N-M  tu 


1»1 


,  b.  +  X 
1  =  1  l 


and  then 


With  the  expression  for  R(X)  in  (7.14)  (or  (7.15)),  repeated  computations  of 
R(X)  for  different  values  of  X  are  straightforward,  once  the  matrix  H  the 
vector  v  and  the  diagonal  of  the  matrix  L  are  computed.  Hence  the  value  of 
X  minimizing  R(X)  can  be  computed  by  standard  minimization  methods. 

We  remark  here  without  proof  that  the  arquments  in  [13)  can  be  used  here  to 
prove  that  the  solutions  to  problems  I-m  and  Il-m  have  representations  of  the 
form 

N  M 

I  c  (x)  ♦  l  d  q  (x) 
j-1  1  1  v»l  V  V 


where 


(7.16) 


</^(x)  =  /  Em(x,t)dt1dt2  ,  t  = 


and  the  {q^}  are  as  before.  The  vectors  c  and  d  satisfy  equations  of  the 
form  (7.4)  and  (7.5)  with  K .  and  T^v  given  by 
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